Correlation effects in the ground state of trapped atomic Bose gases 
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We study the effects of many-body correlations in trapped ultracold atomic Bose gases. We calculate the 
ground state of the gas using a ground-state auxiliary-field quantum Monte Carlo (QMC) method [Phys. Rev. E 
70, 056702 (2004)]. We examine the properties of the gas, such as the energetics, condensate fraction, real-space 
density, and momentum distribution, as a function of the number of particles and the scattering length. We find 
that the mean-field Gross-Pitaevskii (GP) approach gives qualitatively incorrect result of the kinetic energy as 
a function of the scattering length. We present detailed QMC data for the various quantities, and discuss the 
behavior of GP, modified GP, and the Bogoliubov method under a local density approximation. 



I. INTRODUCTION 

The many-body physics in trapped Bose gases has drawn 
intense interest since the experimental realization of Bose- 
Einstein condensation (BEC) in ultracold, dilute alkali 
atoms (TJ] . The systems are "clean" and highly controllable 
experimentally. The dominant interactions are simple and 
well-understood, and the strength of the interatomic inter- 
actions can be readily tuned by means of Feshbach reso- 
nances fai. With the recent realization of degenerate Fermi 
gases jj, 0,01, these ultracold systems provide an ideal "lab- 
oratory" for studying many-body physics. 

In the weakly-interacting regime, mean-field theories work 
quite well, for instance the Gross-Pitaevskii (GP) equation |6i 
0, @1 for boson ground states. Much work has been done 
to study the ground state of the Bose atomic gases beyond 
mean field. For example, a modified GP equation was pro- 
posed 1 9] by inclusion of one-loop quantum corrections and 
the use of local-density approximation. Esry 1 10] developed 
a Hartree-Fock theory as a means of including the correlation 
effects in the BEC many-body calculations. Mazzanti and co- 
workers full applied a correlated basis theory fl2il to study the 
detailed structure of dilute hard- and soft-sphere Bose gases. 
A comparative study for the modified GP and correlated ba- 
sis approaches is presented in Ref. [oj Recently, McKinney 
and co-workers 11411 used a many-body dimensional perturba- 
tion theory to compute the ground-state energy and breathing- 
mode frequency of spherically trapped gases at different inter- 
action strengths. 

Semianalytic methods are versatile and generally very easy 
to extend to realistic systems with large number of particles. 
However, they are approximate and each has its own limita- 
tions, especially in the strongly-interacting regime. Computa- 
tional methods such as quantum Monte Carlo (QMC) provide 
a useful, complementary alternative. A variety of such cal- 
culations have been carried out for atomic boson systems, in- 
cluding variational Monte Carlo fl5ll and the exact diffusion 
Monte Carlo (DMC) (H El 

studies on both the homoge- 
nous 11811 and trapped gases 

We have recently developed an auxiliary-field quantum 
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Monte Carlo (AF QMC) method |22|] for the ground state 
of many-boson systems. While the standard DMC works 
in real space with particle configurations, our method works 
in the second-quantized formalism, which automatically ac- 
counts for particle permutation statistics. The calculation can 
be carried out in any single-particle basis. Conceptually, the 
method provides a way to systematically improve upon mean 
field while retaining its basic machinery, capturing correlation 
effects with a stochastic, coherent ensemble of independent- 
particle solutions. Various observables and correlation func- 
tions can be calculated relatively straightforwardly. 

The initial motivation of this study was to use the AF QMC 
method to quantify, by direct comparison with GP, the effects 
of interactions in trapped Bose gases, and to provide addi- 
tional precise numerical data where they were not available. 
(Although the method is not exact for bosons with repulsive 
interactions, the systematic errors are very small in the pa- 
rameter region of interest, as we show below.) In particular, 
we were interested in the behavior of the system as a function 
of the interaction strength, which, unlike in typical condensed 
matter systems, can be probed directly in experiments. We 
found that GP yielded significant errors in the energetics in the 
Feshbach resonance regime, which resulted in a qualitatively 
incorrect behavior of the kinetic energy in GP as a function 
of the scattering length. To study the origin of these errors, 
we carried out additional calculations using first-order Bogoli- 
ubov results under a local-density approximation (LDA). The 
purpose of this paper is thus to present our QMC data, and 
discuss the behavior of the GP, modified GP, and Bogoliubov- 
LDA methods as benchmarked by QMC. 

The rest of the paper is organized as follows. In Sec. |ll] 
we describe the many-body Hamiltonian. Our QMC method 
is summarized in Sec. [II]] as are the procedures of our GP 
and Bogoliubov-LDA calculations. Results from QMC, GP, 
and first-order Bogoliubov-LDA methods are presented in 
Sec. II VI where we study the energetics of the gas in three di- 
mensions as a function of the number of particle N and the s- 
wave scattering length a s , and examine the density profile and 
momentum distribution. Our study extends to the strongly- 
interacting regime achieveable by Feshbach resonances. In 
Sec. [V] we discuss the implications of our comparisons be- 
tween GP, modified GP, Bogoliubov-LDA and QMC. In addi- 
tion, we also discuss the influence of the details of the two- 
body potential. Concluding remarks are given in Sec. I VII Fi- 
nally, in the appendix, we describe additional details on our 
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Bogoliubov and QMC calculations, including benchmark re- 
sults on the systematic errors in our QMC. 



II. MODIFIED BOSE-HUBBARD MODEL 



We consider N Bose particles in a three-dimensional cube 
of side length 2r b , under the periodic boundary condition. 
Similar to our earlier work 12211 . we use the Bose-Hubbard 
model as the discrete representation of the many-body Hamil- 
tonian on a real-space lattice: 



where 



^=^E fc V(k)^(k) 
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(1) 



where the kinetic energy operator is modified from the Bose- 
Hubbard form we used earlier, and is expressed in momentum 
space instead, with 



0(k) 



1 



(2r fc )3/2 



(2) 



The sum over k is taken over all the (discretized) momentum 
coordinates. Equation Q describes both the homogenous and 
trapped Bose gases. For a homogenous gas, oj q = 0. In both 
cases, we use a large enough r b to minimize the boundary 
effects. We will set h = m = 1 throughout this paper. 

We discretize the cubic simulation box into an L x L x L 
lattice. The lattice spacing is <r = 2r b /L. We enumerate 
the real-space sites using an integral index i ranging from 1 
through L 3 . The coordinate of the i-th site is given by r^. The 
periodic boundary condition restricts the values for the mo- 
mentum coordinates k = (ki, ...,^3) to hi = nrii/rb, where 
m is an integer in the range [— L/2\ < m < [L/2\. We 
will use the index q = 1, 2, L 3 to enumerate the points in 
the momentum space; correspondingly, k g is the momentum 
vector of the q-th point. 

The field operators on the lattice are defined to be 



The discretized Hamiltonian is therefore given by 



b„ 



(3) 
(4) 
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'ho 



(6) 
(7) 



and a ho = yJKfmuj^ is the harmonic oscillator length scale. 
The representation of the kinetic energy in Eq. (0 reproduces 
the continuum spectrum more faithfully than the real-space 
finite-difference form in the original Hubbard form, and al- 
lows quicker convergence with the size of the grid, L. 

The contact two-body potential in the continuum is ill- 
defined ll23l l24ll because of the ultraviolet divergence. The 
momentum-space interaction strength, 



y 2 s(q) = / drV2s(r)e 



-?q-r 



is uniform for any |q|. The discretized Hamiltonian allevi- 
ates the problem to a large degree by introducing a mometum 
space cut-off k c and replacing the (5-potential by an on-site in- 
teraction parameterized by the scattering length, a s . However, 
the discretized two-body potential in Eq. (0 must be adjusted 
in order to yield the correct two-body scattering length, and 
a s in Eq. must be replaced by an appropriate a' s for the 
lattice. Following the standard treatment, we obtain the regu- 
larized a', which for a 3D lattice is ll25ll 



(8) 



1 - 2.442749a s /<r 



For the system to be in the dilute limit and the form of our two- 
body potential to be valid, we need the density at the lattice 
sites to satisfy {hi) <C 1. 



III. COMPUTATIONAL METHODS 
A. Quantum Monte Carlo method 

1. General formalism for many-boson boson ground states 

We briefly describe our method of computing the ground 
state of many bosons. A detailed account can be found in 
Ref.|22[ We project the ground-state wave function |$q) from 
a trial wave function I^O, 



(9) 



where |^ T ) in this study is the GP solution (see Sec. lIiIBl for 
details). The projector 



1 gs — e e 

= e AT&r e^ AT/? e" Arl> e-5 A ^ 



C(Ar 2 



(10) 
(11) 



is evaluated stochastically by rewriting the two-body part into 
a multidimensional integral. 
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The two-body part of the potential in Eq. (0 can be written 
as a sum of the squares of one-body operators V = —^^2ii>f, 
where ij = \f~-U cjc, is essentially the density operator. We 

use the following Gaussian integral identity to rewrite e~~ ArV 
in terms of the one-body operators: 



1 



/2tt 



12 1-2 



e VA^-*)« ; (12) 



where the constant a is determined below. We use an im- 
portance sampling scheme to sample the ground-state wave 
function, so that 



l$o) 



\<P) 



(vf T |0) 



(13) 



where each |</>) is a mean-field solution, i.e., a permanent con- 
sisting of identical single-particle orbitals. In practice, this 
means that each \<p) is represented by a single-particle orbital. 

The projection in Eq. is then realized by random walks 
in a manifold of mean-field solution s l2p.E3l . which are gov- 
erned by the following equation 122112711 : 



\<t>') 



dcrp(er)B{(T - a)W{a, cj>)\<j)) 



where 




TJ e v / A7(Ti-ff I )«i 
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The optimal choice of the constant vector cr is B22LI27 



(* T |0> 



= -VAn 



(14) 



(15) 



(16) 



(17) 



(18) 



With this choice, the weight factor in Eg. dl4> can be written 
in the so-called local energy form I22L 12711 : 



W{(T,(j)) W e - A ^(*T|H|0)/<* T l0) = e -ArB L (0) 



(19) 



In practice, whether the local-energy or the hybrid form in 
Eq. < I17I > is more efficient will depend on the system. For 
the calculations in this paper, we have mostly used the local- 
energy form. 

We initialize a population { \ <p) } to mean-field solutions, 
e.g., |\& T ). A single random-walk step for each walker con- 
sists of updating the orbital and its associated weight w±, 



H°-*)\4>) 



(20a) 
(20b) 



where the auxiliary fields {cr^} are drawn from the Gaussian 
probability density function p(tr). 

The computation of observables is done using the back- 
propagation estimator I22LI26I1 . 
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bp 



(* T | e - T bp H i|$ c 
(* T |e- T bp£|$ ) 



(21) 



which for large enough r bp yields the ground-state expecta- 
tion value for any observable. 



2. Phaseless approximation 

The formalism above is exact. For repulsive interactions, 
unfortunately, ii in Eq. i\2i becomes imaginary. This is sim- 
ilar to the phase problem in fermionic systems l27ll . and we 
apply the recently developed phaseless approximation, which 
has been shown to work well in electronic-structure calcula- 
tions l27Tl . This method eliminates the phase problem at the 
cost of a systematic bias which is dependent on the trial wave 
function. As we will demonstrate in benchmark calculations 
in Appendix|X| the bias is relatively small for the bosonic sys- 
tems we study here. Indeed, for all but the largest values of 
a s , it is possible to perform unconstrained calculations with 
fixed imaginary-time, j3 = nAr, in which j3 can be made 
sufficiently long that essentially exact ground-state values are 
obtained. Comparison with these results shows that the sys- 
tematic error in the phaseless approximation is small (see Ap- 
pendixlATi. 

In the phaseless approximation, the weights {w$} are re- 
stricted to real, positive values. We define the phase rotation 
angle AO by 



A6 = 3 In 



V(* T l</>> 



(22) 



This is the complex-phase rotation of the walker's overlap 
with the trial wave function as a result of the application of 
B(cr — cr) to \(f>). In the phaseless approximation, the evolu- 
tion of w A is altered to 



w. 



cos( A0) I W(<t,< 
0, 



|A0| < it/2 
otherwise 



(23) 



which prevents the walkers from reaching the origin of the 
(^-pl^-complex-plane. Equations J20at and d23l define the 
algorithm of the phaseless QMC method. 

In invoking the phaseless approximation, it is helpful to re- 
arrange the two-body interaction term in H such that a mean- 
field background is subtracted: 

V = ~\ I> - <^» 2 - E W + \ E<^) 2 ' (24) 

i i i 

where the constant (u^) is the mean-field expectation value, 
e.g., with respect to |\P T ). The residual term involving 
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(vi — (vi)) is then used in Eq. H2\ . This would have no ef- 
fect in the exact formalism above, where, as we discussed in 
Ref. I22I the importance sampling transformation effectively 
introduces the background subtraction even if the bare form 
of V is used. With the phaseless approximation, however, the 
rotation angle is controlled by the mixed-estimate of Vi. Re- 
ducing its average by subtracting the mean-field background 
will thus help reduce the rotation, and improve the behavior 
of the approximation in Eq. J23I . 

We note that the presence of phaseless approximation 
breaks the time-reversal symmetry of the ground-state pro- 
jector. The forward, phaseless propagator (e~ Th p H ) p h is no 
longer formally equivalent to the back-propagated, phaseless 

propagator {e~ Tb p H ) ph [see Eq. J21H . This results in an 
additional systematic error in the back-propagation estima- 
tor. The expectation value of an operator A computed from 
back-propagation is (<1>q|A|$q), where |$q) and |$' ') are 
the approximate ground-state wave functions (normalized) in 
the forward- and backward-direction, respectively, and they 
are in general not the same. This is similarly the case in 
the constrained-path Monte Carlo for fermion lattice models 
I2H l28ll . It was shown l28ll that the error vanishes linearly 
as \^ T ) — » l^o)- We will further discuss the effect of the 
phaseless constraint in Sec.[V]and AppendixlAl 



B. GP self-consistent projection and QMC trial wave functions 

We solve the GP equation on the same lattice defined for 
QMC, using a self-consistent projection with the GP propaga- 
tor exp (— AtHqp) l22ll - Aside from a factor (N — 1)/N in 
front of the interaction terms, the one-body Hamiltonian Hgp 
is simply Eq. (0 with the replacement 

c\c\c iCi -> 2(c] Ci ) c\ Ci - (c\ Ci ) 2 , (25) 

where the expectation is with respect to the GP solution. As 
discussed in Ref. f22, our QMC can be thought of as stochas- 
tically carrying out the functional integral, while GP is the 
saddle-point approximation. 

The U parameter in the GP calculations is given by the bare 
a s rather than the regularized a' s using Eq. {HJ, because the 
shape-independent 5 potential has become a mean-field po- 
tential in the GP approximation. It is these GP results that we 
compare with. 

For our QMC calculations, the trial wave function ^ T is 
taken to be the solution of the GP-like projection, but with 
the regularized a' s . This wave function is different from the 
correct GP solution above, which is obtained using the bare 
a s . Each value of the discretization parameter <; corresponds 
to a different renormalized a' s [see Eq. (JHJ], and gives rise to 
distinctly different results, while the correct GP solution con- 
verges rapidly with ^ (see Fig. [8}. As the trial wave function, 
however, we argue that the optimal choice is the best vari- 
ational solution, which is given by the corresponding mean- 
field calculation with the same a' s . 



C. Bogoliubov approximation 

In the Bogoliubov approximation l29l l30l I3TI1 . correlation 
effects are treated by means of perturbation, where the zeroth- 
order term is the GP mean-field solution. The approach was 
first formulated for a homogenous Bose gas. It assumes a 
macroscopic occupancy of the lowest energy state (k = 0), 
namely (N—Nq) <C N. For each density p = N/Cl and inter- 
action strength, the total energy per particle E-gog/N, momen- 
tum distribution 7TB g(k), and condensate fraction Nq/N can 
be written down analytically in the thermodynamic limit. The 
corrections to the mean-field GP results are expressed in terms 
of the gas parameter pa? s , which gives a measure of the devia- 
tion from the mean-field picture. Note that the bare a s should 
be used, since the regularization of the scattering length is im- 
plicitly done in the Bogoliubov approximation as is in GP. 

It is important to truncate the summation over k in com- 
puting the momentum distribution and kinetic energy. This 
stems from the incorrect behavior of the Bogoliubov 7r(k) at 
large momenta: 7r Bog (k) oc l/|k| 4 as |k| — * 00. Physically, 
the form of the two-body potential requires that |k|a s <C 1, 
therefore the contribution from |k| larger than a cutoff mo- 
mentum k c should be excluded. We use an explicit numerical 
summation with the same k-space grid as in QMC. This au- 
tomatically limits the sum to the reciprocal lattice (excluding 
k = 0). In addition, it helps to correlate the finite-size effects 
in the two calculations, and allows for a more direct compari- 
son of the results between Bogoliubov and QMC. 

We extend the Bogoliubov approach to the inhomogeneous 
case using a local-density approximation (LDA), by treating 
each lattice site as a locally homogenous Bose gas. This 
is similar to the LDA approximation for electronic systems 
under density functional theory l32ll . and we refer to it as 
Bogoliubov-LDA. The approximation is expected to be rea- 
sonable if the density is smooth and slowly varying, which is 
fulfilled in our dilute Bose gas systems. 

The kinetic energy, for example, is a sum of two contribu- 
tions under this approach: one from the curvature (inhomo- 
geneity) of the density profile, and the other from Bogoliubov 
correction. Given the real-space density p(r), it is 

(^Bog-LDA =-\j d 3 V^)V 2 ^) 

, (26) 
+ J d 3 rT B og[p(r)]p(r), 

where the functional T-B 0g [p(r)] is the Bogoliubov kinetic en- 
ergy per particle for a gas with uniform density p = p(r). 
More details on our Bogoliubov-LDA procedure are provided 
in AppendixlBl) 



IV. RESULTS 

In this section, we present results on the energetics, con- 
densate fraction, density profile, and momentum distribution. 
Individual energy terms are computed: (T) is the kinetic en- 
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ergy, (V2B) the two-body interaction energy, and (Vtrap) the 
external trapping potential. 

In the calculations, we typically use a 24 x 24 x 24 lattice, 
with a simulation box of linear dimension 2r b = 14a ho . This 
gives us a lattice constant of c = 0.583a ho . Our trap length is 
a ho = 8546 A, which gives typical peak densities of about 10 
to 40 fim~ 3 for 100 to 1000 particles in the trap. The lattice 
constant is large compared to our scattering lengths (up to 
a s ~ 1000 A), which is consistent with the assumption in 
neglecting the details of the two-body potential. 



1 r 

QMC, a s = 80 I 

QMC, a s = 300 ^ 

QMC, a s - 500 h 

GP - 

noninteracting 




FIG. 1: The ground-state column density p y (x,z — 0) of a trapped 
gas containing N = 100 bosons with three different scattering 
lengths: a s = 80 A, 300 A, and 500 A. QMC statistical error bars 
are indicated. The GP densities are next to the corresponding QMC 
curve, and are all shown in dashed lines. Also shown as a reference 
is the non-interacting profile. 



A. Density Profile 

Figure ^ shows the density profiles of 100 trapped bosons 
for three different scattering lengths. To make a connection 
with experiments, we show the column density 



p v (x,z) = / dyp{x,y,z) , 



(27) 



that is, the density integrated along a particular direction (e.g., 
the jy- axis), w hich can be observed through optical measure- 
ments I33II34I35II . As we increase a s , the condensate expands 
due to the increasing repulsive interactions. Similarly, as we 
add more particles into the gas, the density profiles expands, 
as shown in Fig.|2| 

Compared to GP, the QMC peak density is always lowered, 
and the QMC overall density profile is more extended. For 
a s = 80 A, the peak column density is lowered by 0.5% 
from GP. For a s = 500 A, this difference is about 7%. Ear- 
lier many-body calculations using the correlated basis ap- 
proach |ll Hil and DMC (H ifTJ also showed the same 
qualitative behavior. Below we will further discuss these in 
connection with the energetics and momentum distribution. 



t 1 1 1 1 1 

noninteracting, N=50 ■-■ 

GP . 

4 QMC, N = 50 x I 
\ QMC, N = 500 I — * — I 
QMC, N = 1000 o I ' 




x {a J 

FIG. 2: The ground-state column density p y (x,z — 0) of a trapped 
gas of N — 50-1000 bosons with scattering length of a s = 120 A. 




200 400 600 
a, (A) 



200 400 600 

a, (A) 




200 400 600 

a s (A) 



200 400 600 

a s (A) 



FIG. 3: Ground-state energy per particle and its individual compo- 
nents as a function of the scattering length. The system has 100 
bosons. QMC error bars are statistical. 



B. Energetics 

Figure [5] shows the ground-state energy and its individual 
components as a function of the interaction strength. We see 
that, as the scattering length a s is increased, the total energy 
increases as expected. Both GP and Bogoliubov-LDA ener- 
gies are in reasonable agreement with QMC, deviating more 
at larger a s . The GP energy is slightly lower than the exact 
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results (no variational principle due to regularization), while 
Bogoliubov-LDA is higher. The external potential energy, 
(Vtmp), also increases with a s , which is a consequence of the 
expansion of the density profile with interaction, as shown in 
Fig. [0 The GP trap energy is lower than QMC, consistent 
with the result in Fig. that QMC density profiles are more 
extended. 

The kinetic and interaction energies are shown in the bot- 
tom panels of Fig. [3] The discrepancy between GP and QMC 
is more pronounced. In particular, the GP kinetic energy de- 
creases monotonically with a s , because the density profile ex- 
pands and the system becomes less confined. The QMC ki- 
netic energy, on the other hand, shows a nonmonotonic be- 
havior. For small a s , the kinetic energy decreases as a s is in- 
creased, tracking the GP result. At a s > 400 A, however, the 
kinetic energy curves up and increases with a s . The QMC in- 
teraction energy is significantly lower than the mean-field in- 
teraction energy at large a s , and the GP result increases much 
more rapidly with a s than QMC. Indeed the QMC curve ap- 
pears to turn downward at the last point, but our data is not suf- 
ficient to establish this, as it is possible that a larger systematic 
error from the phaseless approximation may have contributed 
to make the QMC result smaller (see the benchmark results in 
Appendix IAV 

From a single-particle picture, we would expect the QMC 
kinetic energy to be lower than that of GP, since the QMC 
density profiles are more extended. In reality, correlation ef- 
fects become more important as a s increases, which raises the 
kinetic energy with interaction. This is illustrated clearly by 
considering the uniform Bose gas, for which we show cor- 
responding results in Fig. |4] The GP ground state is a zero- 
momentum condensate. In the many-body ground state, inter- 
actions excite particles into higher-momentum single-particle 
states, raising the kinetic energy as a result. The QMC results 



physics is qualitatively unchanged if the GP densities are used 
instead. The result shows good agreement with the full many- 
body calculation. In particular, the Bogoliubov kinetic energy 
shows an increase similar to the QMC prediction. The corre- 
sponding interaction energy is also reduced, although not as 
much as in QMC. Overall, the Bogoliubov results capture the 
basic picture and confirm that correlations are an important 
ingredient in the energetics of the gas. 

C. Condensate fraction and momentum distribution 

Figure |5] shows the condensate fraction as a function of in- 
teraction strength. GP by definition gives 100%. We see that 
the actual depletion is about 4% at 800 A. Again, the Bogoli- 
ubov result agrees well with QMC. Figure [6] shows the mo- 




100 200 300 400 500 600 700 800 

a s (A) 



FIG. 5: Condensate fraction as a function of the scattering length. 
The system is the same as that of Fig. [3] The condensate fraction is 
defined as the leading eigenvalue (normalized by N) of the one-body 
real-space density matrix. QMC statistical errors were not shown. 




J l l L 





5. 



100 200 300 400 500 600 

a, (A) 




100 200 300 400 500 600 

a, (A) 



FIG. 4: Kinetic and interaction energies in the uniform Bose gas as 
a function of the scattering length. We show results from QMC, Bo- 
goliubov, and GP. The density is p = 0.542 pm~ 3 . The simulation 
box has 100 particles on a 16 x 16 x 16 lattice, representing a physical 
volume of fi = 184.4 jim 3 . 

in the trapped gas are thus the outcome of the competition be- 
tween mean-field and correlation effects. 

The Bogoliubov-LDA calculations, whose results are also 
shown in Fig. [5] help to quantify this picture further. We 
use QMC density profiles in the calculation (hence the exact 
agreement between the Bogoliubov-LDA and QMC estimates 
of the trap energy in Fig. [3}, although we have verified that the 



mentum distribution for two scattering lengths: a s = 200 A 
and 500 A. The QMC's momentum distribution is more 
peaked than GP. This translates in the real space to a more 
extended density profile for QMC, as is observed in Fig.^ 

The graph also shows the contribution to the kinetic energy 
from various k = |k| regions, since the kinetic energy is re- 
lated to the momentum distribution through 



(f) oc / k 2 dkir(k)k 2 



(28) 



Relative to GP, the QMC distribution is depleted in the 



medium-/s regime, around k ~ a^ . Part of this depletion 
goes to the low-momentum region near k = 0, and the other 
to the high-/c region. At a higher a s , the depletion shifts to- 
ward the smaller k region. It is clear that the enhancement 
in the high-fc region results in the increase of the kinetic en- 
ergy. The kinetic energy is strongly enhanced in the larger a s 
cases, which results in the upturn of the kinetic energy curve 
in Fig. |3] 

A precision measurement of the momentum distribution 
would be useful to reveal the detailed structure of the many- 
body correlations in the Bose gas. Our results from a lattice 
do not have enough resolution to reveal whether there are finer 
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to further improve the MGP equation, and make it more like 
DFT-LDA. For the kinetic energy, however, the MGP would 
give the same qualitative results as GP, even when the ex- 
act xc-functional is used and the exact density is obtained, 
because the "kinetic energy" that is explicitly defined in the 
MGP framework is incomplete. In fact, the same would seem 
to apply to DFT-LDA for electronic systems. This is an im- 
portant conceptual difference between MGP and Bogoliubov- 
LDA approach, although they are closely related and lead to 
the same total energy results. 



FIG. 6: Momentum distribution for trapped gases at two different 
scattering lengths. The system again has N = 100 bosons. The 
left panel shown a cut along the fc-axis. The right panel shows the 
difference between the QMC and GP, multiplied by k 4 . 



structures in the momentum- or real-space density. (A fine 
structure in the density profile was predicted by the DMC cal- 
culations l2lll .1 In the auxiliary-field QMC framework, a bet- 
ter resolution in the density profile may be obtained by choos- 
ing a more suitable basis set, such as Hartree-Fock states fioll . 
whereby the GP solution becomes the lowest-energy state in 
this basis set, and also the leading solution in the ground-state 
wave function. 



V. DISCUSSIONS 



A. GP, modified GP, and Bogoliubov-LDA approaches 



B. Finite-size effects and limitations of the on-site potential 

There are two kinds of finite-size errors in our calculation: 
the error due to finite simulation box size, and the discretiza- 
tion error due to finite lattice constant. The first kind is eas- 
ily reduced, by increasing the simulation box size, ?-{,. In the 
trapped boson calculations with N = 100 particles, we have 
checked that r b > 5a ho is sufficient for a s < 1000. For cal- 
culations with large values of TV, we use r b = 7a ho to al- 
low simulations of large enough condensate while keeping the 
finite-size errors much less than our statistical error. 

The discretization error from the finite lattice constant, <r, is 
more subtle. On the one hand, sufficiently small should be 
used so the results converge to the continuum values. Figure0 
shows the convergence of the total energy. It also illustrates 
the effect of regularizing the scattering length, as discussed 
in Sec.|n] In Fig.|8j we show the convergence of the density 
profile. 



We have shown that the many-body correlations qualita- 
tively change the behavior of the kinetic energy in the trapped 
Bose gas. The Bogoliubov approximation l29l l3(il. l3lll under 
the local density approximation (LDA), which we refer to as 
Bogoliubov-LDA, captures this trend quite well. The LDA 
provides a good way to include the correlation effects based 
on the homogenous Bose gas results. This is perhaps not sur- 
prising, given the diluteness of the gas. 

In contrast, the mean-field GP method by construction ap- 
proximates the kinetic energy only by the part that arises 
from the inhomogeneity of the gas, missing the portion from 
many-body effects. The separation of these two portions is 
especially clear in the homogeneous gas, as we illustrated 
in Sec. IIV Bl This appears to be a rather generic feature 
of independent-particle approaches. The same would ap- 
ply to the modified GP (MGP) method H HHH HI, 
which can be viewed as the bosonic counterpart of the 
standard electronic structure method of LDA under density- 
functional theory (DFT). In that framework, the MGP equa- 
tion is an outcome of using the Bogoliubov results for the uni- 
form Bose gas as the "exchange-correlation" (xc) functional, 
i.e., LDA+Bogoliubov (as opposed to the Bogoliubov-LDA 
above). This method has a great advantage in that it allows 
self-consistent calculations. For example, the real-space den- 
sity can be calculated directly and would not need to be im- 
ported as was done with the Bogoliubov-LDA. Further, it is of 
course possible to use exact QMC results on the uniform gas 
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FIG. 7: The effect of finite discretization on the QMC and GP total 
energies due to the lattice constant <;. The test system has TV = 100, 
a s — 120 A. We show the total energy of the system for <; ranging 
from 0.8 a ho (10 x 10 x 10 lattice) through 0.25 a ho (32 x 32 x 32 
lattice). Also shown is the QMC energy calculated without regular- 
izing the scattering length, which fails to converge. The statistical 
error bars are smaller than the point size. 

On the other hand, the lattice constant is also coupled to the 
on-site potential that we use, which in turn affects the detailed 
energetics of the system. The on-site potential effectively has 
finite range and strength which depend on This is equiv- 
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FIG. 8: The effect of finite discretization on the density profiles in 
QMC and GP. The density profile p(x, y = 0, z — 0) is shown for 
two different values of the lattice constant, The test system has 
TV = 100 and a 3 = 400 A. The GP curves are indistinguishable. 
The densities obtained from the QMC trial wave functions are also 
shown. These are GP-like solutions but with regularized scattering 
lengths. They do not converge like the QMC or true GP densities. 



alent to setting the cutoff momentum k c oc 1/q in the inter- 
action matrix elements. Figure [9] shows the total and kinetic 
energies as <; is varied. The total energy is less sensitive to the 
details of the interaction potential, as are the real-space den- 
sity (see Fig. [8} and the trap energy. The dependence on ? in 
the kinetic and interaction energies, however, is not negligible. 
(This dependence is consistent with the observation of Maz- 
zanti and co-workers lUlll when they varied the range of their 
soft-sphere repulsive potential.) It is important to note that the 
nonmonotonic behavior of the kinetic energy is observed at all 
<; values. As ? is reduced, the upturn is more enhanced, indi- 
cating a stronger effect from the interactions as the potential 
is made narrower and harder. 
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FIG. 9: The effect of the lattice discretization and on-site interaction. 
The total energy is insensitive to <;, whereas the kinetic energy shows 
more dependence on the details of the interaction. The system is the 
same as in Fig. [3] 



Ideally, we would like to decouple the basis-size error (due 
to finite lattice spacing) from the effect of the details of the 
potential. For this purpose, the on-site pseudopotential is in- 
adequate. The (5-function potential is meant to be used with 



the short-distance contributions already "integrated out" |4C 
The effects above represent corrections from the details of the 
interaction potential as defined by the on-site form, which 
change as we vary <r. It is easy to see that in the limit of 
? — ► 0, the gas is trivially noninteracting in the exact many- 
body picture @|, since the range of the interaction potential 
is zero. However, if the conditions for the validity of the po- 
tential are maintained, the corrections should be small and not 
affect essential properties, as we have illustrated. A better 
pseudopotential should have an intrinsic decay in momentum 
space with well-defined convergence properties. 



C. Bias due to phaseless approximation 

The phaseless approximation, as demonstrated by the 
benchmarks in Appendix [X] gives an excellent approxima- 
tion to the true many-body ground state for weak to moderate 
interaction strengths. Nevertheless, systematic errors on the 
computed observables are expected. For example, the varia- 
tional principle, that the total energy computed by QMC is an 
upper bound to the exact energy, is not guaranteed in the pres- 
ence of phaseless approximation 1271. Eill . We even observe 
this bias in the a s = 500 A results shown in Table UTI 

The systematic bias is noticeable, but remains quite small 
up to the largest scattering lengths we study, as can be seen 
from the benchmark data. It is interesting to compare the 
phaseless and unconstrained QMC energies in Table HH At 
a large a s = 500 A, the phaseless approximation lowers the 
kinetic energy (as well as the interaction energy) compared 
to the unconstrained result. This trend is observed for all a s 
values. Since the phaseless bias increases with the interaction 
strength, it should lead to an underestimation of the upturn of 
the kinetic energy. Thus the nonmonotonic behavior of the ki- 
netic energy should actually be slightly stronger than shown 
by QMC. 

We have shown in Ref. |52 that the QMC results is in- 
dependent of the input trial wave function $ T . This is no 
longer the case in the presence of the phaseless approxima- 
tion. The approximation imposes a constraint based on the 
overlap (<3> T |(/)), and each fy T in principle has different con- 
straining properties. This dependence is very weak, however, 
as we observed in our calculations among trial wave functions 
of the same general form (GP-like). 

The phaseless approximation can also affect the Trotter er- 
ror, which arises from the use of a finite time step At in 
Eq. il It . This error is controllable, and can be extrapolated 
away by running at different values of At. Because the ro- 
tation angle in the random walk is proportional to \J At U, 
the severity of the phaseless projection is affected by At, as 
is the extent of the population fluctuation. The latter is im- 
portant in back-propagation, where it is highly desirable to 
keep branching to a minimum. If phaseless projection causes 
significant loss of the population, the Trotter error will be in- 
creased. Procedures that reduce the extent of the phase pro- 
jection, for example, subtracting the mean-field background 
shown in Eq. (1241 . will thus improve computational efficiency 
(in addition to possibly reducing the systematic error). 
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VI. CONCLUSIONS 

We have studied the ground state of realistic systems of 
trapped interacting Bose atomic gases using a many-body 
auxiliary-field QMC method, as well as GP and the Bogoli- 
ubov method under a local density approximation. We ob- 
served the effect of correlations in the energetics, condensate 
fraction, real-space density profiles, and momentum distribu- 
tion. The density profile is more expanded compared to the 
GP prediction. The momentum distribution shows enhance- 
ment in the occupation of the low- and high-momentum states. 
The kinetic energy, contrary to the GP estimate, is not mono- 
tonic with the scattering length a s . The Bogoliubov method is 
able to reproduce this trend qualitatively. Additional calcula- 
tions on the uniform Bose gas were performed to help under- 
stand and quantify our results. 

Through this study we also further tested and developed 
our QMC method. We found that the phaseless approxima- 
tion developed for electronic systems 12711 worked quite well 
in the context of boson calculations with repulsive calcula- 
tions. Because of the simplicity of these bosonic systems 
compared to electronic systems, they have provided an ideal 
testbed and allowed us to carry out more benchmark calcu- 
lations and gain additional insights on controlling the phase 
problem, which is crucial for making QMC more useful for 
a wide variety of problems. It is hoped that the formalism 
we developed will allow the study of many interacting Bose, 
Fermi, and mixed-species systems. The method can also ac- 
count for different external experiment environments ( 1 -D or 
2-D, rotations, anisotropic traps, optical lattices, etc.) quite 
straightforwardly. 
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APPENDIX A: BENCHMARK RESULTS ON THE 
PHASELESS APPROXIMATION IN QMC 

In this appendix, we show benchmark results on the phase- 
less approximation in dealing the complex-phase problem, as 
discussed in Sec. lIII A~2l We first show results on a small sys- 
tem for which exact diagonalization can be done. We choose 
a one-dimensional Bose-Hubbard system. The corresponding 
Gross-Pitaevskii calculation is also done at the same Hubbard 
parameters t, U, and k. (Here U is a fixed parameter which 
is the same in QMC and GP.) Table U compares the energet- 
ics and condensate fraction obtained using various methods: 
exact diagonalization, our QMC with the phaseless approxi- 
mation (ph-QMC), and GP self-consistent projection. 



The ph-QMC improves over GP, and in general agrees well 
with exact diagonalization. The bias due to the phaseless ap- 
proximation is visible in the trap energy (Vtrap)- In our phase- 
less QMC calculation, the mean-field background was sub- 
tracted in the Hamiltonian, as shown in Eq. J24I . Applying 
the phaseless approximation directly leads to more population 
fluctuations in the random walk and larger systematic errors 
in (Vtrap) and (V 2 b)- 

TABLE I: Benchmark of QMC with the phaseless approximation 
(ph-QMC) against exact diagonalization. The test system is a ID 
Bose-Hubbard Hamiltonian with 13 sites and 4 particles. The param- 
eters are t = 2.676, U = 1.538, and « = 0.3503. QMC statistical 
errors are in the last digit, and are shown in parantheses. Error bar in 
the condensate fraction was not estimated. 

Type E (f) (V trap ) (V 2B ) (N /N) 

Exact 4.244 1.183 1.793 1.268 98.5% 

ph-QMC 4.242(8) 1.182(6) 1.799(1) 1.262(3) 98.4% 

GP 4.429 1.029 1.800 1.599 100.0% 



We now show calculations on a large system with realistic 
a s values. We use the unconstrained QMC (u-QMC) as the 
reference. For weak to moderate interaction strength, the un- 
constrained QMC can be carried out for a short period of time 
r before the signal is completely lost in large Monte Carlo 
fluctuations. To obtain the desired accuracy, we perform many 
short QMC runs and average the results. For each scattering 
lengths, we verified that the short runs have reached conver- 
gence with respect to the projection time. The severity of the 
phase problem grows rapidly with a s , and such runs are not 
possible for large values of a s . 

Table|n]shows the phaseless QMC with the local-energy ap- 
proximation [Eq. J19H for 3D trapped gas of 100 atoms with 
a s — 80 A and 500 A. The first case represents a typical situ- 
ation in the trapped atomic gas experiments far from Feshbach 
resonances, while the second is a medium-strength interaction 
deep into the range of a s we study. The At parameter was ad- 
justed so that the Trotter error is similar to or smaller than the 
statistical error. We see that the agreement between the phase- 
less and unconstrained calculations is good. 

As a further check, we compare our QMC result on the uni- 
form Bose gas with an earlier diffusion Monte Carlo (DMC) 
calculation by Giorgini and co-workers fljil . which is exact. 
We use their results for the soft sphere potential with large ra- 
dius of R = 5a s , which best matches our situation, namely 
<T ~ 2R ~ 10a s . As we show in the left panel of Fig.^| our 
results agree well with their DMC energies. 



APPENDIX B: BOGOLIUBOV GROUND STATE 

The Bogoliubov approximation for the homogenous Bose 
gas assumes a macroscopic occupancy of the lowest energy 
state (k = 0), namely (N — N ) < N. We will work in 
the thermodynamic limit, N — > oo and f2 — > oo, keeping 
the density p = N/Cl finite. The creation and annihilation 
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TABLE II: Benchmark of QMC calculations with and without the 
phaseless constraint for a s = 80 A and 500 A. We simulate 100 
atoms in a 3D harmonic trap with a ho — 8546 A. The simulation 
lattice is 24 x 24 x 24. Shown here are per-particle quantities. All 
energies are expressed in the unit of hco . 

Type E/N (f)/N {Vt^ P )/N (V2b)/N~ 

a s = 80 A 

ph-QMC 1.7943(3) 0.5984(3) 0.96029(9) 0.23562(8) 
u-QMC 1.7947(2) 0.5987(2) 0.96006(4) 0.23594(4) 

GP 1.7924 0.5947 0.95649 0.24121 

a 3 = 500 A 

ph-QMC 2.6777(2) 0.500(3) 1.5638(6) 0.591(1) 
u-QMC 2.6811(4) 0.511(7) 1.563(2) 0.614(3) 
GP 2.620 0.408 1.4901 0.721 



The summation, however, must be performed with care, as 
mentioned in Sec. HIT J The analytic results for the energy 
and condensate fraction, Eqs. JB4> and JB9L are obtained by 
extending the summation variable to infinity, because the con- 
tribution from outside the ka s <C 1 region is assumed to be 
small. This assumption does not hold for the kinetic energy, 
since the sum diverges due to the unphysical nature of n(k) at 
large |k|. 

To benchmark our Bogoliubov approach, we perform QMC 
and Bogoliubov calculations for a homogenous Bose gas at 
different scattering lengths, as shown in Fig. ^3 We com- 
pute the energetics and condensate fraction using three differ- 
ent methods: GP, Bogoliubov, and QMC. As we see here, the 
first-order Bogoliubov approximation estimates the energet- 
ics and condensate fraction very well for a small enough gas 
parameter (here pa^ < 10 -4 ). 



operators for the zero-momentum state are approximated as 
scalars, 



0t(o)«0(O) 



(Bl) 



We then ignore all terms higher than quadratic in the remain- 
ing creation and annihilation operators. The form of the two- 
body potential also requires that ka s <C 1. Within this ap- 
proximation, the energy per particle is given by l4lll 



-Efiog = E Bog /N 

Airpa, 
= 2N 



lirpa s 



1 



150F 

and the occupation of the k momentum state by |42 




™Bog( k ) = 



l-al 



(MO). 



where 
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(B3) 

(B4) 
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The quantity £ = (Snpa s ) 1 I 2 is the healing length lEoll of 
the condensate. The condensate fraction is given by 



^ 1 1 n,\ 

3V 7T ' 

The kinetic energy per particle can be computed through 



N ~ N 
= 1 - 



(B8) 
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FIG. 10: Benchmark of our QMC and Bogoliubov-LDA calculations 
in the uniform Bose gas. The system is the same as that in Fig.|4| 
The upper triangle data points in the g.s. energy plot are from diffu- 
sion Monte Carlo (DMC) calculations using a soft sphere (SS) poten- 
tial fiUl . In the condensate fraction plot, we also show the analytic 
Bogoliubov result without truncation of |k|, Eq. <B9> . 



Note that the condensate fraction estimated by Bogoliubov 
with the truncation in the sum in fc-space agrees much better 
with QMC than the analytic Bogoliubov. The analytic Bogoli- 
ubov estimate is off, as discussed above, because it is extrapo- 
lated to an infinite box size, and it includes contributions from 
very high momentum states. 

We note that the kinetic energy, which is very small in the 
small a s regime, is no longer negligible for larger a s values. 
For a s = 600 A, or equivalently paj? = 1.2 x 10 -4 , the kinetic 
energy (see Fig. |4ji is about 37% of the total energy. This 
is consistent with our discussion in Sec. I1VI on the balance 
between the mean-field and correlation effects. 

We can extend the Bogoliubov analysis above to deal with 
the case of a inhomogeneous, trapped gas. We use the so- 
called local-density approximation (LDA) by treating each 
lattice site as a locally homogenous gas. The density profile 
p(r) can be estimated using GP or any other methods which 
provides a good approximation to the density profile. Using 
the same fc-space lattice as QMC, we compute the "local" en- 
ergetics (per particle) and condensate fraction. The density 
is then used to weight-average the local contributions. The 
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Bogoliubov-LDA estimate of the kinetic energy is 

COeog-LDA =~\f d 3 r^) V 2 ^p~(v) 

/ (Bll) 

The interaction energy is given by 

(^2B>Bog-LDA = J d 3 r (£ Bog [p(r)]-f[p(r)] Bog )p(r). 

(B12) 



[1] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, 
and E. A. Cornell, Science 269, 198 (1995). 

[2] S. L. Cornish, N. R. Claussen, J. L. Roberts, E. A. Cornell, and 
C. E. Wieman, Phys. Rev. Lett. 85, 1795 (2000). 

[3] B. DeMarco and D. S. Jin, Science 285, 1703 (1999). 

[4] F. Schreck, G. Ferrari, K. L. Corwin, J. Cubizolles, 
L. Khaykovich, M.-O. Mewes, and C. Salomon, Phys. Rev. A 
64, 011402 (2001). 

[5] A. G. Truscott, K. E. Strecker, W. I. McAlexander, G. B. Par- 
tridge, and R. G. Hulet, Science 291, 2570 (2001). 

[6] E. P. Gross, Nuovo Cimento 20, 454 (1961). 

[7] E. P. Gross, J. Math. Phys. 4, 195 (1963). 

[8] L. P. Pitaevskii, Sov. Phys.-JETP 13, 451 (1961). 

[9] E. Braaten and A. Nieto, Phys. Rev. B 56, 14745 (1997). 
[10] B. D. Esry, Phys. Rev. A 55, 1147 (1997). 
[11] F. Mazzanti, A. Polls, and A. Fabrocini, Phys. Rev. A 67, 
063615 (2003). 

[12] S. Fantoni and A. Fabrocini, in Microscopic Quantum Many- 
Body Theories and Their Applications, edited by J. Navarro 
and A. Polls (Springer- Verlag, Berlin, 1998), Lecture Notes in 
Physics Vol. 510, p. 119. 

[13] A. Fabrocini and A. Polls, Phys. Rev. A 60, 2319 (1999). 

[14] B. A. McKinney, M. Dunn, and D. K. Watson, Phys. Rev. A 69, 
053611 (2004). 

[15] J. L. DuBois and H. R. Glyde, Phys. Rev. A 63, 023602 (2001). 
[16] D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980). 
[17] C. J. Umrigar, M. P. Nightingale, and K. J. Runge, J. Chem. 

Phys. 99, 2865 (1993). 
[18] S. Giorgini, J. Boronat, and J. Casulleras, Phys. Rev. A 60, 5129 

(1999). 

[19] D. Blume and C. H. Greene, Phys. Rev. A 63, 063601 (2001). 
[20] A. R. Sakhel, J. L. DuBois, and H. R. Glyde, Phys. Rev. A 66, 
063610 (2002). 

[21] J. L. DuBois and H. R. Glyde, Phys. Rev. A 68, 033602 (2003). 



The external trap energy is straightforward to compute, 
namely 



(Wbos-lda = / d 3 ry trap (r>(r). (B13) 



[22] W. Purwanto and S. Zhang, Phys. Rev. E 70, 056702 (2004). 
[23] N. P. Proukakis, K. Burnett, and H. T. C. Stoof, Phys. Rev. A 

57, 1230 (1998). 
[24] B. D. Esry and C. H. Greene, Phys. Rev. A 60, 1451 (1999). 
[25] Y. Castin, J. Phys. IV France 116, 89 (2004), also available at 

http : //arXiv . org/abs/cond-mat/0407118 
[26] S. Zhang, J. Carlson, and J. E. Gubernatis, Phys. Rev. B 55, 

7464 (1997). 

[27] S. Zhang and H. Krakauer, Phys. Rev. Lett. 90, 136401 (2003). 
[28] J. Carlson, J. E. Gubernatis, G. Ortiz, and S. Zhang, Phys. Rev. 

B 59, 12788 (1999). 
[29] N. N. Bogoliubov, J. Phys. (U.S.S.R.) 11, 23 (1947). 
[30] T. D. Lee, K. Huang, and C. N. Yang, Phys. Rev. 106, 1135 

(1957). 

[31] T. T. Wu, Phys. Rev. 115, 1390 (1959). 

[32] W. Kohn, Rev. Mod. Phys. 71, 1253 (1999). 

[33] M. R. Andrews, M.-O. Mewes, N. J. van Druten, D. S. Durfee, 

D. M. Kurn, and W. Ketterle, Science 273, 84 (1996). 
[34] M. R. Andrews, D. M. Kurn, H.-J. Miesner, D. S. Durfee, C. G. 

Townsend, S. Inouye, and W. Ketterle, Phys. Rev. Lett. 79, 553 

(1997). 

[35] L. V. Hau, B. D. Busch, C. Liu, Z. Dutton, M. M. Burns, and 

J. A. Golovchenko, Phys. Rev. A 58, R54 (1998). 
[36] A. Fabrocini and A. Polls, Phys. Rev. A 64, 63610 (2001). 
[37] G. S. Nunes, J. Phys. B: At. Mol. Opt. Phys. 32, 4293 (1999). 
[38] A. Banerjee and M. P. Singh, Phys. Rev. A 64, 63604 (2001). 
[39] H. Fu, Y. Wang, and B. Gao, Phys. Rev. A 67, 053612 (2003). 
[40] A. J. Leggett, Rev. Mod. Phys. 73, 307 (2001). 
[41] K. Huang, Statistical Mechanics (John Wiley & Sons, New 

York, 1987), 2nd ed. 
[42] The continuum momentum distribution, which is often referred 

to in the main text, per particle, is given by 7r(k) /N = w (k) = 

n(k)/[(2 7 r) 3 (0 ]. 



